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Abstract 

Let a measurement consist of a linear combination of damped complex exponential modes, plus 
noise. The problem is to estimate the parameters of these modes, as in line spectrum estimation, vibration 
analysis, speech processing, system identification, and direction of arrival estimation. Our results differ 
from standard results of modal analysis to the extent that we consider sparse and co-prime samplings 
in space, or equivalently sparse and co-prime samplings in time. Our main result is a characterization 
of the orthogonal subspace. This is the subspace that is orthogonal to the signal subspace spanned 
by the columns of the generalized Vandermonde matrix of modes in sparse or co-prime arrays. This 
characterization is derived in a form that allows us to adapt modern methods of linear prediction and 
approximate least squares, such as iterative quadratic maximum likelihood (IQML), for estimating mode 
parameters. Several numerical examples are presented to demonstrate the validity of the proposed modal 
estimation methods, and to compare the fidelity of modal estimation with sparse and co-prime arrays, 
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versus SNR. Our calculations of Cramer-Rao bounds allow us to analyze the loss in performance sustained 
by sparse and co-prime arrays that are compressions of uniform linear arrays. 

Index Terms 

Co-pime array, IQML, modal analysis, orthogonal subspaces, sparse array 

I. Introduction 

In this paper, we investigate the problem of estimating the parameters of damped complex exponentials 
from the observation of non-uniform samples of their weighted sum. This problem arises in many 
applications such as modal analysis, speech processing, system identification, and direction of arrival 
(DOA) estimation. 

There is a vast literature on different modal estimation methods from uniformly sampled time or space 
series data, starting with the work of Prony 0. Other methods include approximate least squares or 
maximum likelihood estimation 0, 0, reduced rank linear prediction 0, 0, MUSIC 0, and ESPRIT 
1171 . 0. While there are extensions of MUSIC and ESPRIT for direction of arrival estimation from 
non-uniformly sampled data (see, e.g., ll9l-| fT3l ), Prony-like methods have mainly been developed for 
uniformly sampled data, and extending such methods to non-uniformly sampled data has not received 
much attention (exceptions being lfl4l and JT3J). 

Non-uniform sensor array geometries, without aliasing ambiguities, have a long history in sensor array 
processing, dating back to minimum-redundancy arrays lfl6l . The introduction of co-prime arrays in lfT71 . 
m, and llT9l has created renewed interest in such geometries. In this paper, we consider two specific 
cases of non-uniform sensor arrays. These are sparse arrays and co-prime arrays. Both of these geometries 
can be viewed as subsampled (or compressed) versions of a dense uniform line array, whose consecutive 
elements are separated by a half wavelength in space[j] Specifically, the sparse array can be thought 
of as a subsampled version of a dense uniform line array, plus an extra sensor that is positioned at a 
location on the array that allows us to resolve aliasing ambiguities. The co-prime array consists of two 
uniform subarrays, each obtained by uniformly subsampling a dense uniform line array with co-prime 
subsampling factors. The co-prime property allows for resolving aliasing ambiguities. 

'if we were sampling in time, then the dense sequence of uniform samples would have had spacings equal to the Nyquist 
interval. 
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Naturally, any subsampling in space results in a reduction in signal-to-noise ratio (SNR), by the com¬ 
pression factor, and leads to a loss in estimation performance. Our studies in BOl . ll2Tll . and Il22l address 
the effect of compression on Fisher information, the Cramer-Rao bound, and the probability of a swap 
between signal and noise subspaces. Assuming that the loss in SNR due to compression has tolerable 
effects on estimation or detection, or can be compensated by collecting more temporal snapshots (requiring 
a scene to remain stationary for a longer period), the question is how can methods of linear prediction 
and approximate least squares be adapted to the estimation of mode parameters in sparse and co-prime 
arrays? In this paper, we address this question. 

We determine a parameterization of the orthogonal subspace. This is the subspace that is orthogonal 
to the signal subspace spanned by the columns of a generalized Vandermonde matrix of the modes in 
sparse and co-prime arrays. This parameterization is of a form that is particularly suitable for utilizing 
approximate least squares, such as iterative quadratic maximum likelihood (IQML) (see |[2j], (3|, and 
fl23l ). for estimating the modes. Although we present our numerical results in the context of sensor array 
processing, all of our results apply to the estimation of complex exponential modes from time series 
data. Our numerical results here, and in ll20l . ll2H . and ll22l , show that there is a loss in performance 
sustained by sparse and co-prime arrays that arc compressions of uniform linear arrays. A rough rule of 
thumb is that effective SNR is reduced by 101og 10 C', where C is the compression ratio. For example, 
in our experiments a 50-element array is subsampled to a 14-element co-prime array, for a compression 
ratio of 50/14. The loss in SNR is roughly 5.5 dB. 

Remark 1: A small number of other authors have also considered estimating the parameters of complex 
exponentials from non-uniformly sampled data using Prony-like methods. In iTBl . the authors approach 
the modal estimation problem by fitting a polynomial to the non-uniform samples and estimating the 
parameters of the exponentials using linear regression. For the case that the modes are on the unit circle, 
in na a truncated window function is fitted to the non-uniform measurements in the least squares sense, 
and then an approximate Prony method is proposed to estimate the frequencies of the exponentials. 
These approaches are different from ours and do not involve characterization of orthogonal subspaces 
for utilizing modern methods of linear prediction. 
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II. Problem Statement 


Consider a non-uniform line array of m sensors at locations I = {io, i\, ■ ■ ., i m -i} i n units of half 
wavelength in space. We assume, without loss of generality, that iq = 0. Suppose the array is observing 
a weighted superposition of p damped complex exponentials (modes). These modes arc determined by 
the mode parameters Z}. = Pkt J<>k , k = 1,2,... ,p, where the /;:th mode has a damping factor pi- and an 
electrical angle 0}- £ (—7r, 7r]. Suppose the array collects N temporal snapshots. Then, the measurement 

equation for the Ith sensor (located at ifi can be written as 

p 

yi[n\ = y^ x k[n\z l k + e;[n], n = 0,1,..., N - 1, (1) 

fc=i 

where n is the snapshot index, :c/.[n] denotes the amplitude (or weight) of the A th mode at index n, and 
ei[n] is the measurement noise at sensor l. In vector form, we have y[n] £ C m , 


y[n] = V(z, I)x[n] +e[n], n = 0,1,... ,N - 1, 


( 2 ) 


where y[n] = [yo[n],yi[n\,... ,y m -i[n\\ r is the array measurement vector, x[n] = [xi[n] : x' 2 [ri ], 
...,x p [n]] T is the vector of mode amplitudes at index n, e[n] = [eo[n], ei [n],..., e m _i[n]] T is the 
noise vector at index n, and V(z,I) £ C mxp is a generalized Vandermonde matrix of the modes 
z = [z 1} z 2 , ■ ■ ■, z p ] T , given by 


V (z, I) = 


po 


po 


y*0 


(3) 




We consider the case where x[n] is free to change with n, and assume that the e/ [n] ’s, are i.i.d. complex 
normal with mean zero and variance cr 2 . This means that the measurement vectors y[n], n = 0,1,..., iV— 
1 are i.i.d proper complex normal with mean V(z,I)x[n] and covariance cr 2 l. Under this measurement 
model, the least squares estimation and the maximum likelihood estimation of the modes {zk] p k= i and 
mode weights {x[/t] are equivalent and can be posed as 


mm 
z,x[0].x[AI—1] 


AT-1 

£ 

71=0 


y[n] - V(z,I)x[n] 


(4) 


The least squares estimate of x[n] is 


x[n] = V + (z,I)y[n], 


(5) 
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where V + (z,I) = (V' ff (z,I)V(z,I)) 1 V i/ (z,I) is the Moore-Penrose pseudoinverse of V(z,I). The 

least squares estimate of the modes is obtained as 

JV-l 

z = argmin ^ y ff [n](I - P V (z,i))yN 

n=0 

N-l 

= argnun ^ y H [n]P A ( Z)I) y [n ], (6) 

n =0 

where A(z,I) is a full column rank matrix that satisfies 

A^(z, I)V (z, I) = 0 (m _ p)xp , (7) 

and Pv( z ,i) an d Pa(z,i) = I—Pv(z,i) are the orthogonal projections onto the column spans of V(z, I) and 
A(z,I), respectively. We denote these column spans by the subspaces (V(z,I)) and (A(z,I)). We call 
(V(z,I)) the signal subspace and (A(z,I)) the orthogonal subspace. Note that (A(z,I)) = (V(z,I))- L . 
See Figure [T] 



Fig. 1. The signal subspace (V(z,I)) and the orthogonal subspace (A(z,I)) = (V(z,I)) X - In the figure, we have dropped 
(z, ,1) and have simply used A, V, (A), and (V). 


For a given array geometry, the basis matrix V(z, I) given in ([3]), and the subspace (V(z, II)), are fully 
characterized by the p modes z = [z\, Z 2 , ■ ■ ■, z p ] T . This subspace, parameterized by z, is an element of 
a Grassmanian manifold of dimension p. Now, let us rewrite V(z,I), using elementary operations, and 
with some abuse of notation, as 


V(z,I) 


Vr(z,I) 

V 2 (z,I) 


(8) 
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where Vi(z,I) E C pxp is invertible and V 2 (z,I) E C 1 ”"' p) x p. Then the basis matrix A(z,I) for the 
orthogonal subspace is the Hermitian transpose of 

A h (z,I) = [-V 2 (z,I)V 1 - 1 (z,I) | I m _ p ], (9) 

Although these p-dimensional characterizations of the signal and orthogonal subspaces have minimum 
parameterization z E C p , it is not easy to solve the least squares problem ([6]) using these characterizations. 

For an m-element uniform line array, a particular /^-parameter characterization of A(z, I) exists that makes 
solving ([6]) relatively simple |[T|. We will review this characterization in Section [111] Then, we derive such 
suitable parameterizations of A(z, I) for two specific non-uniform arrays: sparse and co-prime. 

• Sparse array. In this case, the location set I is given by I s = {0, d, 2 d, ..., (m — 2 )d,M}, where M 
and d are co-prime integers, that is, (M. d) = 1, and d > 1. This array may be thought of as two 
subarrays. The first is a downsampled version, by a factor d, of an (m — 1)(/-element uniform line 
array (ULA) with half wavelength interelement spacings. The second is a single sensor at location 
M in the line array such that M and d are co-prime. We call this the sparse array because of the 
single element that sits apart from the origin of the first subarray. We note that M need not be 
greater than (m — 2 )d. 

• Co-prime array: In this case, I = Ii U I 2 , where Ii = {0, m 2 , 2m 2 ,..., (mi — l)m 2 }, I 2 = 
{mi, 2mi,..., (2m 2 — l)mi}, and (mi,m 2 ) = 1. Again the array is composed of two subarrays. 
The first is an mi-element ULA with interelement spacings of m 2 and sensor locations Ii. The 
second is a (2m 2 — l)-element ULA with interelement spacings of mi and sensor locations I 2 . This 
co-prime geometry was recently introduced in ffTTIl and lfl8l . 

Remark 2: In both cases, the co-prime constraint guarantees that aliasing ambiguities due to undersampling 
can be resolved. Although a sparse array can be viewed as a special case of a co-prime array, we consider 
them separately, because it is easier to first derive a suitable characterization of the orthogonal subspace 
(A) for a sparse array, and then generalize it to a co-prime array. Our parameterizations are not minimal. 
They involve 2 p parameters, instead of p, but as we will show in Section [TV] they are specifically designed 
to utilize modern methods of linear prediction and approximate least squares, such as IQML. 

Remark 3: By now it should be clear that A(z,I), and therefore its parameterization, depend on both 
the mode vector z and the array geometry I. Therefore, from here on, we may drop (z,I) and simply 
use A, V, (A), and (V). 
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III. Characterization of the Orthogonal Subspace for Uniform Line Arrays 


Consider a uniform line array of m equidistant sensors located at l u = {0,1, 2,..., m — 1}, taking 
measurements from the superposition of p modes as in ([2]). The signal subspace in this case is characterized 
by the Vandermonde matrix V in <|3]) with I = I u . To characterize the orthogonal subspace (A), consider 
the polynomial A(z)\ 


A(z) = - ZkZ -1 ) 

k =l 

v 

= '^a i z~ l \ a 0 = 1 (10) 

i =o 

which has (zi,Z 2 ,, z p ) as its p complex roots. The (m — p) dimensional orthogonal subspace (A) is 


spanned by the m — p linearly independent columns of A: 

CLp CLp —1 ■ ■ ■ oq 1 0 


A = 


H 


ijp ,p- 

0 a r 


jjp ... ... Ur 1 


( 11 ) 


l 0 • • • 0 a v 

Since A ;/ V = 0, and the columns of A are linearly independent, V and A span orthogonal subspaces 
(V) and (A) in C m . The above parameterization is at the heart of methods of linear prediction, 
approximate least squares, and IQML (see, e.g., lf2ll—lf5ll). 

Using the p-parameter representation for (A) in (JTT|) we may re-write the least squares problem of ^ 


as 


N—l 

a = argmin E y"HPAyH. 

a=[a 1 ,...,a p ] T eCP n=0 


( 12 ) 


There are many algorithms to approximately solve the nonlinear least squares problem in (12). One 
approach is to ignore the (A a A) _1 term in the projection matrix Pa = A(A /; A)~ 1 A ;/ and solve the 


following modified least squares or linear prediction problem: 

AT-l 

a = argmin y^[n]AA tf y 
aeo ^ 


n 


(13) 


n =0 


The iterative quadratic maximum likelihood (IQML) algorithm (see O, 0, and f23l ) is another method 
to approximately solve ( [T2j ). In the /th iteration of IQML, the parameters a i are estimated by iteratively 
minimizing the quadratic form 


a i = argmin a l 

a; £C P 


H 


'N -1 


^2 YH [ n \( A i 1 -iA l -i)- 1 Y[n\ 


n =0 


a i, 


(14) 
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where A/_i is formed as in ( [IT] ) using the estimated parameters a;_i from iteration (l — 1) and y[n] is 
the following (m — p) x p Hankel data matrix for snapshot n: 

yo[n\ ■■■ y P -i[n\ y p [n\ 

yi N ••• y P [n] y P + tN 


Y\n] = 


J/m-l-pM 


ym- 2 [n] y m -i[n] 


5 U'P J 


(15) 


. The 


After a number of these iterations the sequence {a/} converges to an estimate a = [di,. 
polynomial A(z) = Y^ =(J a t z^ 1 is then formed from this estimate and its roots are taken as the mode 
estimates (z\, Z 2 , ■ ■ ■, z p ). 

IV. Characterization of the Orthogonal Subspaces for Sparse and Co-prime Arrays 

In this section, we present simple characterizations of the orthogonal subspace (A) for the sparse and 
co-prime arrays discussed in Section [IIJ Based on our characterizations, we adopt IQML for approximate 
least squares estimation of complex exponential modes in such arrays. 

A. Sparse Array 


Consider the sparse array described in Section |TTJ The set of sensor locations for this array is Jl, s = 
{0, d, 2d ,..., (m — 2 )d, M}. The generalized Vandermonde matrix V in this case is 


V(z,I s ) = 


(m—2)d (m—2)d 


n 


Zo 


Z■ 


(m—2)d 


P 


~M 

Z 1 


z 2 




(16) 


For d > 1 it is clear that without the use of the last sensor at location M, we cannot unambiguously 
estimate the modes, because any two modes Zk and Zke^ 2nq ^ d , q = 1 , 2 , ■ ■ ■ ,d — 1 produce the same 
measurement. This is the aliasing problem for subsampled arrays. 

To characterize the (m — p) -dimensional orthogonal subspace (A), determined by the modes {zkYi ,, 
we first form the polynomial A(z) from the r/th powers of Zk, namely the = zf„ k = 1,2,... ,p: 


A{z) = \\{l - w k z 1 ) = J2 


diZ 


a 0 = 1. 


(17) 


k =1 


i =0 
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Since {wifc}f =1 are the roots of A(z), the first m — p — 1 columns Ao of A £ p), w hj c h j s to 

satisfy A^V = 0 , can be written as 


A h - 
— 


dp dp —i 

0 On 


Ol 1 0 


JLp - - • ... n t f 0 


(18) 


0 • • • 0 a v 

But of course any mode of the form z k e^ 2 ’ Kq ^ d , q = 1..... d—1, would produce the same w k and therefore 
the same Aq. This is the ambiguity caused by aliasing. 


Now, consider the polynomial 


B(z) = z M + Y,biZ^- i)d . 


(19) 


i=l 


Suppose the coefficient vector b = [ 6 i, 62 , ■ ■ ■ , b p ] T is such that the actual modes {z k } p k=l arc the roots 
of B(z). That is, B(zk) = 0 for k = 1, 2,... ,p. Then, since M and d are co-prime, for 1 < q < d — 1 
and 1 < k < p we have 


B{z k e j2nq/d ) = z M e i^ Mq /d + J2 b iZ 

i =1 

= z M( e ^ Mq /d - 1) 


(p—i)d 

k 


/ 0 for q = 1 , 2 ,..., d — 1 . 


( 20 ) 


Therefore, the only common roots of B(z), and the dth roots of {w k } p =1 , are {z k Y k =v which are the 
actual modes to be estimated. In this way, B(z) resolves the ambiguities. 


Now suppose {'(Cfc }? =1 are known (or estimated). Then from (191, b can be found by solving the linear 
system of equations 

/ 1 f 1 \ 


b p b p -1 ... bi 


Vi 


„(P~ l) d 


Z n 




= _ 1 ~M M 
\ Zi z 2 




( 21 ) 


which if zf / z d for i / j (as we assume) has a unique solution. 
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Using the 2 p coefficients {a,} and { b, Y‘_ ( , we can characterize (A) by writing A G C mx ("' p) as 

i H 

a p G’n —i " * ■ up 1 0 ■ ■ ■ 0 


A = 


Lp Up - 

0 a r . 


0 • • • 0 a 

bp b p —i • • • b\ 0 


ip - • • up 1 0 


0 1 


( 22 ) 


To estimate a = [up, ■ ■ ■ , a p ] T and b = [ftp, • • • , b p ] T , we need to solve the following problem: 


N-l 


nnn 

a,b 


Y y H N p A y W, 


(23) 


n=0 


We approximate the solution to this problem in two steps. First, we ignore the last column of A and 
estimate a as 


N-l 


a = argmin Y y^N p A r ,yN- 


(24) 


n =0 


In the noiseless case, where the y[n], n = 0,1,..., N — 1, lie in (V), it can be shown that if m > 


2p + 1 then the solution to (24 1 is unique and yields the coefficients of the polynomial A(z) with roots 
(wi,W 2 , ■ ■ ■, w p ). See Appendix A. 


The minimization problem in (24) can be solved using IQML. Now, given a, we form the polynomial 

p 

A(z) = 1 + Y^ biiZ~ l 


i —1 


k= 1 


W k Z 


-U 


(25) 


and derive its roots as {w k } p k=1 - But we know from the structure of the problem that w k = z^, and any 
of the d-th roots of z'l is a candidate solution. Therefore, we construct the candidate set, 

n = {(z ie ^/ d , z 2 e j2wq2/d , • • • , z v e^ qp ' d ) |0 < q u q 2 ,. .., q p < d - 1}, (26) 

which contains all modes and their aliased versions. 

In the second step, to find the p actual modes and resolve aliasing ambiguities, we solve the following 
constrained linear prediction problem: 

N-l 

b = argmin Y Iz/m-iN + C T u[n]| 2 

^ n =0 


s.t. B^( z) = 0, z G 72., 


(27) 
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where u[n] = \yo[n], y\ [n],... ,y p _i[n]] T , and the polynomial B^(z) is obtained from replacing b by £ 
in (Q9J). 


In the noiseless case, where y[n], n = 0,1,.. 
and yields the actual modes. See Appendix B. 


, ./V — 1 lie in (V), the solution b to (21) satisfies (21) 


Our algorithm for estimating modes in a sparse array may be summarized in the following steps: 


1) Estimate a = [fii, & 2 , • • •, a p ] T from (24) using IQML; 

2) Root A(z) to return roots {wk} p k=l - Then, recognizing that the dth roots of vy. are z\ t e^‘ 1 ' Kq l d for 


some q £ {0,1, 2,..., d — 1}, form the set of candidate modes 1Z as in (26); 


3) Solve (27 1 for b; 


4) Intersect the roots of B(z) with 72. 


B. Co-prime Array 


Consider an rn = m ± + 2 7772 — 1 element co-prime array, consisting of two uniform subarrays: one with 
mi elements at locations Ii = {0, 7712 , 2 m 2 ,..., ( 777-1 ~ 1 ) 777 - 2 } and the other with 2 rri 2 — 1 elements 
at locations I 2 = {mi,2mi,..., (2 to 2 — l)mi}, where ( 7711 , 7772 ) = 1 and 7771 > 7772 . In this case, the 
generalized Vandermonde matrix V £ C mxp of modes may be partitioned as 


V(z,Ii) 
V(z,I 2 ) 


(28) 


where 



1 

1 

1 



Z 1 

„m 2 

y rn 2 

^p 


V(z,Ii) = 

2m 2 

Z 1 

2m 2 

z 2 

Jm 2 

^ p 

(29) 


(m 1 — l)m 2 
Z 1 

(m 1 — l)m 2 

Z 2 

(m 1 — l)m 2 

• Zp 
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and 

.mi mi mi 

z1 *2 *p 

~2 m, 2)nx 2mi 

V (z, I 2 ) = _ 2 _ P (30) 

(2m 2 — l)m 1 (2m 2 —l))77i (2m 2 — ljtri! 

Z 1 ^2 ' ' ' Z p 

are the Vandermonde matrices for the two individual subarrays of the co-prime array. 

Let Ai € C miX ( mi-p ) and Bi € C( 2m2_1 ) x ( 2m2_1_p ) be matrices that are orthogonal to V(z,Ii) and 
V(z, I 2 ), respectively. That is, V(z, Ii) = 0 and V(z, I 2 ) = 0. Following our results in the sparse 
case, we may parameterize Ai £ C m ' x( ™' 1 ~p) as 

d p (Ip —| ... CL\ 1 ... Q 

0 a p 0 

0 • • • 0 a p ■■ ■ d\ 1 

where {di\ p i=l are the coefficients of a polynomial A(z), whose roots arc Wk = z™ 2 , k = 1 . 2 ,..../). 
That is, 

p 

A(z) = J|(l - Wkz -1 ) 
k=\ 
p 

= '^2 a iZ~ l , do = 1. (32) 

i=0 

Similarly, we parameterize Bi £ (£(2m 2 -i)x(2m 2 -i-p) as 

bp bp —1 • • • b\ 1 • • • 0 

0 b p ■ ■ ■ • • • 0 

0 • • • 0 b p ■■ ■ b\ 1 

where {b l } 1 ‘ =x arc the coefficients of a polynomial B(z), whose roots arc = zT 1 , k = 1.2,..../). 
That is, 

p 

b ( z ) = na _ Skz ~ l ) 

k=\ 

p 

= Y J biZ~\ 6o = l. (34) 

i =0 
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Note that we still need p more independent columns to fully characterize the basis matrix A for the 
orthogonal subspace (A). However, using our partial characterization, we can estimate the modes (with 
no aliasing ambiguities) in the following steps: 

1) Separate the measurements of the two subarrays as u[n] = {y, [n ]|i E Hi} and v[n] = {y.;[n] \i E I 2 }; 

2) Estimate a = [di, a 2 ,..., a p ] T using IQML on u[n]; 

3) Root A(z) to return the roots {■«;/,.}} \- Then, recognizing that the 771 , 2 th roots of up- are Zkf- ]2lTq ' m2 
for some q E {0,1,, m 2 — 1}, form the set of candidate modes 72 .1 as 

Tlx = {(zxe j2nkl/m2 ,z 2 e j2nk2 / m2 ,- ■ • , z p e j27rk ^ m2 ) | 0 < k u k 2 ,..., k p < m 2 - 1}; (35) 

4) Estimate b = [bx,b 2 ,... ,b p ] 2 using IQML on v[n]; 

5) Root B(z) to return the roots {sfc}? =1 . Then, recognizing that the 777 , 1 th roots of .§/, are %e j27r9//mi 
for some q E {0,1,..., m\ — 1}, form the set of candidate modes 1Z 2 as 

72 2 = {(zxe j2 ^ m \z 2 e j2nk2/mi ,- • • , z p e j27rkp ^) \ 0 < ki,k 2 , • •., k p < mi - 1}; (36) 


6) Intersect TZx and Tl 2 , in other words look for the closest (based on the Euclidean metric) p members 
of the set 72 .1 to the set Tl 2 . 

Remark 4: To complete the 2p— parameterization of the basis matrix A for the orthogonal subspace (A), 
consider the standard representation of A given in Define Vi = V(z,I Pl ) and Vo = V(z,I P2 ) where 
I p , = {0, 7712 ,..., (p — 1 ) 771 - 2 } and I P2 = {r?7i, 2mi,... ,pmx}. Then, from ([9]) the p remaining columns 
of A may be represented in Ci E C nxp as: 


C'l [^0 I ^px(m 1 —p) I Ip I Opx(2?7i 2 — 1— p)] 

where 0/,. x /. denotes a k x l matrix with zero entries, I p is the p x p identity matrix, and 

= -V(z,Ip 2 )V~ 1 (z,Ip 1 ) E C pxp . 


(37) 


(38) 


From (37 1 and (38 1 we can see that C\ only depends on {2™'}} j and (z™ 2 }} , which arc obtained from 


a and b by rooting A a (z) and Bb(z) in (32 1 and (34 1 , respectively. Therefore, the full, and minimally 


parameterized, characterization of the orthogonal subspace for the co-prime array may be written as 
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A = 


Ai 


0 


Cl 


(39) 


0 B l 

We note that we do not need this full characterization for estimating the modes. The partial characterization 
using Ai and Bi suffices, at the expense of p fitting equations. 


V. Numerical Results 

In this section we present numerical results for the estimation of damped complex exponential modes 
in co-prime, sparse and uniform line arrays. We consider a ULA of 50 elements. We form our co-prime 
and sparse arrays with 14 elements by subsampling this ULA. For the sparse array, we subsample the 
measurements of the ULA by a factor of d = 4 and place a sensor at M = 3. For the co-prime array, the 
first subarray includes mi = 7 elements with interelement spacing of m 2 = 4, and the second subarray 
includes 2m-2 — 1 = 7 elements with interelement spacing of mi = 7. 

It is insightful to first look at the beampatterns of sparse, co-prime, and uniform line arrays for the 
problem of estimating undamped modes. In this case, the beam pattern B(6) is 

m —1 

B{d) = J2 eji ‘ e ■ (4°) 

z=o 

Figure [2] shows the beam patterns for different array geometries. Although the co-prime and sparse arrays 
of 14 elements have the same aperture and the same main lobe width as the ULA with 50 elements, we 
see that they suffer from higher sidelobes, suggesting that there will be performance losses in resolving 
closely spaced modes using these arrays, relative to the ULA. 

Let us also look at numerical results for the Cramer-Rao bound (CRB) associated with the co-prime, 
sparse and uniform line arrays (See Appendix C). Figure [3] shows the CRB in the estimation of the mode 
z\ = 1 in the presence of an interfering mode z-z- The per sensor SNR is 10 dB. As the interfering mode 
zz gets closer to z\, the CRB in the estimation of the z\ increases. The CRB for the sparse and co-prime 
arrays are similar, but they are higher than the CRB for the ULA. Since the aperture of the three arrays 
are equal, the fewer number of sensors in the sparse and co-prime arrays can be considered the only 
reason for this difference in the CRB. 

We now consider the performance of the approximate least squares estimation methods, shown in Figs. [4]- 
[H The two modes to be estimated here are zi = e j052 and z 2 = 0.95e ?a69 . We choose the per sensor 
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Fig. 2. Beampatterns for ULAs with 14 and 50 elements, a sparse array with 14 elements, d = 4, and M = 3, and a co-prime 
array with 14 elements, m i = 7, and m 2 = 4. 


SNR values for the sparse and co-prime arrays to be 5 dB higher than the SNR for the ULA, based 
on our insight in ll22l about the threshold SNR for ULA and co-prime arrays. When SNR is decreased, 
there comes a point where some of the components of the orthogonal subspace better approximate 
the measurements than some of the components of the signal subspace. This leads to a performance 
breakdown in the estimation of the modes. The SNR at which this catastrophic breakdown occurs is 
called the threshold SNR (see If24l and G2lD . For the compression ratio of 50/14, the threshold SNR for 
the co-prime and sparse arrays is almost 5 dB more than its value for the ULA, which is a consequence of 
the subsamplings by these arrays. We emphasize that any compression increases the SNR threshold. The 
use of a co-prime or sparse array instead of a dense uniform line array is only justified in applications 
where SNR is high enough for the desired estimation resolution, or when SNR can be built up from 
temporal snapshots using long observation periods. The latter of course requires the scene to remain 
stationary over the longer estimation period. 
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CM 

N 


N 




(a) ULA 



0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 



(b) Sparse array 



0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 



(c) Co-prime array 

Fig. 3. The CRB in dB for estimating z\ = 1 in the presence of an interfering mode z 2 : (aj ULA with 50 elements, (bj sparse 
array with 14 elements d = 4 and M = 3. (c) co-prime array with 23 elements mi = 7 and m 2 = 4. For all arrays per sensor 
SNR is 10 dB. 
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(a) (b) 

Fig. 4. Estimating two closely spaced modes z\ = e-' 0,52 and zi = 0.95e J °' 69 using a ULA with 50 elements: (a) Per sensor 
SNR = 0 dB. (b) Per sensor SNR = -5 dB. 




(a) (b) 

Fig. 5. Estimating two closely spaced modes z\ = e 2 °' 52 and Z 2 = O.OSe- 70 ' 69 using a sparse array with 14 elements, d = 4 
and M = 3: (a) Per sensor SNR = 5 dB (b) Per sensor SNR = 0 dB. 
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180 




♦ Actual Modes 

♦ Aliased Modes: First Subarray 
□ Aliased Modes: Second Subarray 

60 


120 


150 


150 


210 


330 


210 


330 


240 


270 


270 


300 


♦ Actual Modes 

♦ Aliased Modes: First Subarray 
□ Aliased Modes: Second Subarray 

60 


(a) 


(b) 


Fig. 6. Estimating two closely spaced modes z i = e-' 0 ' 52 and 22 = 0.95e- f0 ' 69 using a co-prime array with 14 elements, 
mi = 7 and m 2 = 4: (a) Per sensor SNR = 5 dB (b) Per sensor SNR = 0 dB. 


Appendix A 


Let m > 2 p and zf 7 ^ zj for i ^ j. Also, let a = (ai, 02, ..., a p ) J be the solution to (24 1 . That is, 


iv -1 


a = argmin ^ y"[n]P Ao (a)yN, 


(41) 


n =0 


where Ao(cc) denotes an Ao of the form (18 1 , with a% s replacing a*’s for i = 1, 2,... ,p. We wish to 
show that in the noiseless case the roots of the polynomial A(z) = 1 + Yfd=i a i z ~ l are w k = zf for 
k = 1,2,... ,p (the dth power of the actual modes). 


Without loss of generality we assume that N = 1. Let y = y [0]. In the noiseless case, the minimum 
value of the objective function y ff PA 0 ( Q )y i s zero, and because Aq (a.)Ao(a) is always full rank, for 
the solution vector a, we have A^(a)y = 0. Now, in the noiseless case, y = V(z,I s )x, where V(z, I s ) 
is given by ( fT6| ) and x = \x\, X 2 , ■ ■ ■, x p ] T is the vector of mode weights. Therefore, we have 

Ajf (a) V (z, I s )x = 0, (42) 
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which we can reorder to get 

1 1 


~d 

Z 1 


~d 

z 2 


(m—p—l)d (m—p—l)d 


(m—p—l)d 


X\ 0 

0 X 2 

0 0 



A{w 1 ) 


A(w 2 ) 


A(w p ) 


= 0. 


(43) 


Because the matrix on the left hand side of ( |43| ) is a full column rank Vandermonde matrix, and the 
diagonal matrix in the middle is nonsingular (by the assumption of having p actual modes), the above 
equality holds iff 

A (wk) = 0 for k = 1,p. (44) 

Appendix B 


Let m > 2p and zf ^ zf for i / j. Also, let /3 be the solution to (27) in the noiseless case. Here we 


use (3 instead of b for the solution to distinguish noiseless and noisy cases. We show (3 solves (21 1 and 
therefore resolves aliasing. 

Again, without loss of generality, we assume N = 1. Based on our argument in Appendix A, in the 
noiseless case we have 


n = {( Zl e j2nkl/d , z 2 e j2nk2ld , • • • , z p e j2lTk r /d ) \ 0 < k h k 2 ,..., k p < d - 1 }, 


where {zi} p i=l are the actual modes. In this case (27) can be rewritten as: 


(3 = argmin \y m -i + C u| s.t. V p (rj)C = ~r) , rj £ 1Z, 

where rj = ( 771 , rj 2 ,... , r/ p ), r) &M = [r/f 7 , ..., rj^] T and 


(45) 


(46) 


Vpfa) = 


vf 


vt 


rip 


(47) 


Jp-W Jp~ l ) d ... Ap- 1 ) 11 

'll 'h 'Ip 

In the noiseless case, u = V p (z)x, y m -1 = 1 x i z i* ar *4 V p (r 7 ) = V p (z). Therefore, we have 

p 


\y m -i + f3 T u\ 2 = | ^ Xizf 1 + /3 r V p (z)x| 
i =1 

= 1 ]r>( 2 "-■,")!> 0. 


(48) 


i= 1 
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Now, because (3 is the solution to (461, then in the noiseless case \y m -i + /3 1 u | 2 = 0 and from (481 we 
have -q' )M = z® M almost surely. Therefore, 


(3 = —(v£(77))-V m 
= —(Vp(z)) _1 z 0M 

= b, (49) 


where b is the solution to (|2T|). ■ 


Appendix C 


The Fisher information matrix for the measurement model in ([2]) may be written as 


N 


j(z) = y;j B (z), n = 0,1,... ,N — 1, 


n= 1 


where J n (z) is the Fisher information matrix for the estimation of the modes z = [zi, z-i, 
y[n] in ([2]). That is, 


J n (z) = ^G^(z)G n (z), 


o 


where G n (z) = [gi [ra], g 2 [n], • • • ,g p [n}}, and 


Sim = 


i 0 z t 


*0 — 1 


A z i 


Ji —l 




Xi ra 


(50) 


, z p ] T from 


(51) 


(52) 


is the Zth sensitivity vector, for !</</>. The CRB for the estimation of the A th mode zy- is the A:th 
diagonal element of .1 1 (z). 
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